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LONG-TERM  GOALS 

The  long-term  technical  goal  is  to  design,  develop  and  test  the  next  generation,  primitive  equation, 
ocean  model  for  high-resolution  scientific  (ROMS:  Regional  Ocean  Modeling  System)  and  operational 
(TOMS:  Terrain- following  Ocean  Modeling  System)  applications.  This  project  will  improve  the  ocean 
modeling  capabilities  of  the  U.S.  Navy  for  relocatable,  coastal,  coupled  atmosphere-ocean  forecasting 
applications.  It  will  also  benefit  the  ocean  modeling  community  at  large  by  providing  the  current 
state-of-the-art  knowledge  in  physics,  numerical  schemes,  and  computational  technology. 

OBJECTIVES 

The  main  objective  is  to  produce  a  tested  expert  ocean  modeling  framework  for  scientific  and 
operational  applications  over  a  wide  range  of  spatial  (coastal  to  basin)  and  temporal  (days  to  seasons) 
scales.  The  primary  focus  is  to  implement  the  most  robust  set  of  options  and  algorithms  for  relocatable 
coastal  forecasting  systems  nested  within  basin-scale  operational  models  for  the  Navy.  The  system 
includes  some  of  the  analysis  and  prediction  tools  that  are  available  in  Numerical  Weather  Prediction 
(NWP),  like  4-dimensional  variational  data  assimilation  (4D-Var),  ensemble  prediction,  adaptive 
sampling,  and  circulation  stability  and  sensitivity  analysis. 

APPROACH 

The  structure  of  TOMS  is  based  on  ROMS  because  of  its  accurate  and  efficient  numerical  algorithms, 
tangent  linear  and  adjoint  models,  variational  data  assimilation,  modular  coding  and  explicit  parallel 
structure  conformal  to  modem  computer  architectures  (both  cache-coherent  shared-memory  and 
distributed  cluster  technologies).  Currently,  both  ROMS  and  TOMS  are  identical  and  continue 
improving  and  evolving.  ROMS  remains  as  the  scientific  community  model  while  TOMS  becomes  the 
operational  community  model. 

ROMS/TOMS  is  a  three-dimensional,  free-surface,  terrain-following  ocean  model  that  solves  the 
Reynolds-averaged  Navier-Stokes  equations  using  the  hydrostatic  vertical  momentum  balance  and 
Boussinesq  approximation  (Haidvogel  et  al.  2000,  2008;  Shchepetkin  and  McWilliams,  2005,  2009). 
The  governing  dynamical  equations  are  discretized  on  a  vertical  coordinate  that  depend  on  the  local 
water  depth.  The  horizontal  coordinates  are  orthogonal  and  curvilinear  allowing  Cartesian,  spherical, 
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and  polar  spatial  discretization  on  an  Arakawa  C-grid.  Its  dynamical  kernel  includes  accurate  and 
efficient  algorithms  for  time-stepping,  advection,  pressure  gradient  (Shchepetkin  and  McWilliams 
2003,  2005),  several  subgridscale  parameterizations  (Durski  et  al.,  2004;  Warner  et  al.,  2005)  to 
represent  small-scale  turbulent  processes  at  the  dissipation  level,  and  various  bottom  boundary  layer 
formulations  to  detennine  the  stress  exerted  on  the  flow  by  the  bottom. 

Several  adjoint-based  algorithms  exist  to  explore  the  factors  that  limit  the  predictability  of  the 
circulation  in  regional  applications  for  a  variety  of  dynamical  regimes  (Moore  et  al.,  2004,  2009). 
These  algorithms  use  the  ideas  of  Generalized  Stability  Theory  (GST)  in  order  to  identify  the  most 
unstable  directions  of  state-space  in  which  errors  and  uncertainties  are  likely  to  grow.  The  resulting 
singular  vectors  can  be  used  to  construct  ensembles  of  forecasts  by  perturbing  initial  and  boundary 
conditions  (optimal  perturbations)  and/or  perturbing  surface  forcing  (stochastic  optimals).  Perturbing 
the  system  along  the  most  unstable  directions  to  the  state-space  yields  information  about  the  first 
(ensemble  mean)  and  second  (ensemble  spread)  moments  of  the  probability  density  function.  Given  an 
appropriate  forecast  skill  measure,  the  circulation  is  predictable  if  low  spread  and  unpredictable  if 
large  spread. 

There  are  several  adjoint-based  variational  data  assimilation  algorithms  available  in  ROMS.  For  cases 
in  which  the  dynamics  are  imposed  as  a  strong  constraint  ( i.e .  no  model  error  assumed),  there  is  an 
incremental  4DVar  (Moore  et  al.,  2010a;  Powell  et  al.,  2008)  driver  (I4D-Var)  similar  to  that  used 
operationally  at  some  numerical  weather  prediction  centers.  In  the  case  where  errors  are  admitted  in 
the  model,  there  is  an  indirect  representer-based  weak  constraint  4DVar  (Di  Lorenzo  et  al.,  2007; 
Muccino  et  al.,  2008)  driver  (R4D-Var)  and  a  weak  constraint  Physical  Space  Analysis  System  (4D- 
PSAS)  driver. 

There  are  several  biogeochemical  models  available  in  ROMS.  In  order  of  increasing  ecological 
complexity  these  include  three  NPZD-type  models  (Franks  et  al.,  1986;  Powell  et  al.,  2006;  Fiechter  et 
al.,  2009),  a  nitrogen-based  ecosystem  model  (Fennel  et  al.,  2006,  2008),  a  Nemuro-type  lower  level 
ecosystem  model  (Kishi  et  al.,  2007),  and  a  bio-optical  model  (Bissett  et  al.,  1999). 

ROMS  includes  a  sediment-transport  model  with  an  unlimited  number  of  user-defined  cohesive  (mud) 
and  non-cohesive  (sand)  sediment  classes  (Warner  et  al.,  2008).  Each  class  has  attributes  of  grain 
diameter,  density,  settling  velocity,  critical  stress  threshold  for  erosion,  and  erodibility  constant.  A 
multi-level  bed  framework  tracks  the  distribution  of  every  size  class  in  each  layer  and  stores  bulk 
properties  including  layer  thickness,  porosity,  and  mass,  allowing  the  computation  of  bed  morphology 
and  stratigraphy.  Also  tracked  are  bed-surface  properties  like  active-layer  thickness,  ripple  geometry, 
and  bed  roughness.  Bedload  transport  is  calculated  for  mobile  sediment  classes  in  the  top  layer. 

ROMS  is  a  very  modem  and  modular  code  written  on  F90/F95.  It  uses  C-preprocessing  to  activate  the 
various  physical  and  numerical  options.  The  parallel  framework  is  coarse-grained  with  both  shared- 
memory  (OpenMP)  and  distributed-memory  (MPI)  paradigms  coexisting  in  the  same  code.  Because  of 
its  construction,  the  parallelization  of  the  adjoint  is  only  available  for  MPI.  Several  coding  standards 
have  been  established  to  facilitate  model  readability,  maintenance,  and  portability.  All  the  state  model 
variables  are  dynamically  allocated  and  passed  as  arguments  to  the  computational  routines  via 
dereferenced  pointer  structures.  All  private  arrays  are  automatic;  their  size  is  determined  when  the 
procedure  is  entered.  This  code  structure  facilitates  computations  over  nested  grids.  There  are  three 
types  of  nesting  capabilities  in  ROMS:  (i)  refinement  grids  which  provide  increased  resolution  (3:1  or 
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5:1)  in  a  specific  region;  (ii)  mosaics  which  connect  several  grids  along  their  edges,  and  (iii)  composite 
grids  which  allow  overlap  regions  of  aligned  and  non-aligned  grids. 

WORK  COMPLETED 

One  of  our  major  accomplishments  during  the  last  year  was  the  overhaul  of  the  4D-Var  data 
assimilation  capabilities.  ROMS/TOMS  uniquely  supports  three  different  4D-Var  methodologies: 
incremental  strong  constraint  4D-Var  (I4D-Var),  a  physical-space  statistical  analysis  system  (4D- 
PSAS),  and  a  representer-based  4D-Var  (R4D-Var).  To  our  knowledge,  ROMS/TOMS  is  the  only 
open-source,  ocean  community  modeling  framework  supporting  all  these  variational  data  assimilation 
methods  and  other  sophisticated  adjoint-based  algorithms. 

The  vertical  terrain-following  coordinate  in  ROMS/TOMS  was  further  generalized  to  support  an 
additional  transformation  and  multiple  vertical  stretching  functions.  This  development  is  described  in 
https://www.myroms.org/wiki/index.php/Vcrtical_S-coordinatc  in  the  ROMS  wiki.  The  vertical 
stretching  function,  C(ct),  is  one-dimensional,  nonlinear,  continuous,  and  monotonic.  It  is  defined  in 
terms  of  a  few  input  parameters  and  the  local  grid  bathymetry.  This  allows  the  user  to  increase  the 
vertical  resolution  where  it  is  most  needed  in  a  particular  application.  In  addition,  the  new  vertical 
transfonnation  reduces  the  amount  of  numerical  diapycnal  diffusion  in  the  horizontal  advection  of 
tracers. 

We  held  two  very  successful  workshops.  The  first  was  held  at  the  Jean  Kuntzmann  laboratory,  Saint 
Martin  d'Heres  Campus,  Grenoble,  France  (Oct  6-8,  2008;  http://www.myroms.org/2008workshop). 
The  second  was  held  at  the  Sydney  Institute  of  Marine  Science,  Sydney,  Australia  (Mar  3 1  -  Apr  2, 
2009;  http://www.myroms.org/2009workshop).  As  in  the  past,  several  tutorials  were  offered  in  each 
workshop  about  basic  and  advanced  ROMS  algorithms.  The  next  workshop  will  be  held  at  the  Hawaii 
Imin  International  Conference  Center,  University  of  Hawaii  at  Manoa,  Honolulu,  Hawaii,  April  5-8, 
2010  (http://www.myroms.org/2010workshop).  We  will  also  organize  a  Data  Assimilation  Workshop 
at  the  University  of  California  Santa  Cruz,  Santa  Cruz,  CA,  summer  2010.  This  will  be  a  one  week 
workshop  that  will  include  lectures  on  4D-Var  data  assimilation  methodology  and  algorithms. 

RESULTS 

The  4D-PSAS  and  R4D-Var  algorithms  can  be  configured  as  strong  constraint  (dynamical  model  is 
error  free)  or  weak  constraint  (errors  are  admitted  in  dynamical  model).  All  these  algorithms  are  fully 
working,  parallel  and  affordable.  If  the  same  a  priori  hypothesis  about  the  errors  is  assumed  in  a  strong 
constraint  application  (Moore  et  ai,  2010  a,  b),  they  produce  identical  solutions  and  the  computational 
cost  required  to  identify  the  solution  is  the  same  in  each  case.  This  is  discussed  in  detail  by  Courtier 
(1997)  and  can  be  used  to  check  if  the  ROMS/TOMS  4D-Var  algorithms  are  coded  correctly.  Figure  1, 
shows  the  4D-Var  minimization  cost/penalty  function  J,  a  measure  of  the  misfit  between  model  and 
observations,  for  a  4-day  assimilation  cycle  (1-4  July,  2000)  in  a  regional  model  of  the  California 
Current  System  (described  in  Broquet  et  a/.,  2009).  The  left  panel,  Fig.  la,  shows  the  value  of  J  for 
I4D-Var  (red  curve,  model  space  search)  and  strong  constraint  4D-PSAS  (blue  curve,  observation 
space  search).  After  35  inner  loop  iterations,  both  algorithms  converge  to  the  same  J  value.  The  rate  of 
convergence  is  also  the  same.  The  theoretical  minimum  value  Jmm  of  J,  as  described  in  Courtier 
(1997)  dual  formulation,  is  also  plotted  (black  dashed-curve).  The  fact  that  J  does  not  asymptote  to 
this  theoretical  minimum  indicates  that  one  or  more  of  our  a  priori  error  hypothesis  was  not  entirely 
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correct.  Recall  that  our  prior  hypothesis  about  the  model  error  is  that  the  model  is  perfect  which  is 
unlikely  to  be  true.  The  right  panel,  Fig.  lb,  shows  similar  values  for  J  but  for  the  representer-based 
R4D-Var  algorithm  configured  as  strong  constraint  (blue  curve)  and  weak  constraint  (black  solid 
curve).  The  values  for  I4D-Var  (red  curve)  are  repeated  here  for  comparison.  As  in  Fig.  la,  both 
algorithms  converge  to  the  same  value  and  rate  after  35  iterations.  However,  the  weak  constraint 
solution  (black  solid  curve),  which  accounts  for  errors  in  the  dynamical  model,  is  much  closer  to  the 
Jmin  predicted  by  the  theory.  Indeed,  we  know  that  dynamical  models  actually  have  errors.  These 
errors  tend  to  grow  and  affect  predictability  of  the  ocean  circulation  (forecast  skill). 

xlQ4  xIO4 


Iteratio  Iteration 


Figure  1:  4D-Var  minimization  cost/penalty  function,  J,  versus  number  of  inner  loop  iterations  for 
a  4-day  assimilation  cycle  in  a  regional  model  of  the  California  Current  System  (1-4  July,  2000):  (a) 
comparison  between  strong  constraint  I4D-Var  (red  curve)  and  4D-PSAS  (blue  curve),  (b) 
comparison  between  strong  constraint  I4D-Var  (red  curve),  and  strong  (blue  curve)  and  weak  (black 
solid  curve)  constraint  R4D-Var.  In  all  cases,  the  initial  conditions,  surface  forcing,  and  boundary 
conditions  are  adjusted.  The  dashed  line  indicates  the  theoretical  minimum  value  Jmm  of  J. 

This  is  a  very  important  result  indicating  that  the  4D-Var  algorithms  in  ROMS/TOMS  are  correctly 
coded.  We  get  identical  solutions  regardless  of  weather  the  minimization  is  done  in  model  space  search 
(I4D-Var)  or  observation  space  search  (4D-PSAS  and  R4D-Var).  These  identical  results  are  also 
shown  in  Figure  2.  The  panels  show  maps  of  sea  surface  temperature  (SST)  increments  5xsst(0),  to  the 
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(a)  I4D-Var  (b  4D-PSAS  (c)  R4D-Var 

Figure  2:  4D-Var  sea  surface  temperature  (°  C)  initial  conditions  increments,  Sxsst(0),  to  background  or  first  guess  value,  xb  sst(0), 
from  previous  assimilation  cycle  for:  (a)  strong  constraint  I4D-Var,  model  space  search,  (b)  strong  constraint  4D-PSAS,  observation 
space  search,  and  (c)  strong  constraint  R4D-Var.  The  sea  surface  temperature  increment  patterns  are  identical  for  all  three  4D-Var 

data  assimilation  methodologies. 
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background  initial  conditions  ,  xbsst(0),  for  (a)  strong  constraint  I4D-Var,  (b)  strong  constraint  4D- 
PSAS,  and  (c)  strong  constraint  R4D-PSAS  .  Recall  that  the  new  analysis  is  xa(0)=xb(0)+  5x(0),  where 
xb(0)  is  the  background  initial  conditions  and  5x(0)  is  the  4D-Var  assimilation  increment. 

The  maps  in  Figure  2  are  identical  within  round  off.  Notice  that  the  data  assimilation  corrected  the  sea 
surface  temperature  associated  with  upwelling  along  the  California  Coast. 

All  critical  components  common  to  each  4DVar  approach  including  the  conjugate  gradient, 
preconditioning,  and  error  covariance  modeling  were  greatly  improved.  Coding  and  testing  these 
algorithms  was  technically  challenging  and  time  consuming.  The  conjugate  gradient  algorithm  was 
completely  rewritten  for  all  three  4D-Var  algorithms.  It  is  now  based  on  the  Lanczos  method,  as 
described  in  Fisher  (1998).  Two  versions  of  Limited  Memory  Preconditioners  (LMP)  are  available. 
These  are  the  spectral  preconditioner  as  described  in  Fisher  (1998),  and  the  Ritz  preconditioner  as 
described  in  Tshimanga  et  al.  (2008).  The  technical  description  of  these  algorithms  and  how  to  use 
them  is  documented  in  https://www.myroms.org/wiki/index.php/Variational_Data_Assimilation  in  the 
ROMS  wiki.  The  error  covariance  modeling  was  enhanced  by  imposing  a  multivariate  constraint  via 
the  balance  operator,  K,  as  reported  in  Weaver  et  al.  (2005).  The  state  vector  is  split  between  balanced 
and  unbalanced  components,  except  for  temperature,  which  is  used  to  establish  the  balanced  part  of  the 
other  state  variables.  This  allows  extracting  information  about  the  unobserved  variables  from  the 
observed  data.  It  models  the  off-diagonal  terms  of  the  error  covariance,  B  =  KB„KT,  where  Bu  is  the 
unbalanced  error  covariance  matrix.  Additionally,  new  algorithms  were  developed  to  estimate  the 
posterior  error  or  confidence  in  the  resulting  circulation  estimates  (Moore  et  al.,  2010  a,b). 

IMPACT/APPLICATIONS 

This  project  will  provide  the  ocean  modeling  community  with  a  freely  accessible,  well-documented, 
open-source,  terrain-following,  ocean  model  for  regional  nowcasting  and  forecasting  that  includes 
advanced  data  assimilation,  ensemble  prediction,  and  analysis  tools  for  adaptive  sampling  and 
circulation  dynamics,  stability,  and  sensitivity. 

TRANSITIONS 

The  full  transition  of  ROMS/TOMS  to  the  operational  community  is  likely  to  occur  in  the  future. 
However,  the  ROMS/TOMS  algorithms  are  now  available  to  the  developers  and  scientific  and 
operational  communities  through  the  website  https://www.myroms.org. 

RELATED  PROJECTS 

The  work  reported  here  is  related  to  other  already  funded  ONR  projects  using  ROMS.  In  particular, 
the  PI  (H.  Arango)  closely  collaborates  with  A.  Moore  (adjoint-based  algorithms)  at  University  of 
California,  Santa  Cruz,  A.  Miller  and  B.  Cornuelle  (ROMS  adjoint  and  variational  data  assimilation)  at 
Scripps  Institute  of  Oceanography,  E.  Di  Lorenzo  (Southern  California  predictability)  at  Georgia 
Institute  of  Oceanography,  and  J.  Wilkin  (Mid- Atlantic  Bight  variational  data  assimilation)  at  Rutgers 
University. 

The  PI  is  also  supported  by  the  following  grants: 
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“Characterization  and  Modeling  of  the  Philippine  Archipelago  Dynamics  Using  the  ROMS  4DVAR 
Data  Assimilation  System”,  grant  number  N00014-04-1-0417,  https://www.myroms.org/philex. 
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